In [233]:
from IPython.display import Image
from IPython.core.display import HTML 
In [234]:
Image(filename= "headline.png")
Out[234]:

Introduction

  • Antimicrobial resistance is becoming a worldwide problem

  • Many pathogens are becoming resistant for many different types of antibiotics

  • The institute of advanced studies (IAS) is doing research to the spreading of gonorrhoea in Amsterdam

  • For gonorrhoeae there might be a new antibiotic available within years

  • It is thus important to find a suitable management plan to prevent new resistance

Mathematical model

  • Aim of our research was to build a network model to model the disease spreading
  • Our network models are based on an existing mathematical model by Fingerhurth et al. (2006) [1]
    • Simple network model
    • Stochastic network

Why use a network model?

  • Sexual network portrays the sexual inter-relationships within a group of people and it is possible to study the spread of STIs in the context of social interactions [2]
  • ODEs, in contrast, fail to capture the complexity of connections between cases
  • By using a network a heterogeneous population can be introduced, which adds more detail in the properties of individuals. This can for example be useful in identifying high-risk groups.
  • Adding this heterogeneity can have influence on the disease dynamics [3]
  • Network properties can for example derived from data, in Amsterdam this is for example possible with GGD data

Mathematical model of antibiotic resistance N. gonorrhoeae spread

Below an image is shown which summarized the ODEs of [1] model in a scheme. Shown are the three different classes, susceptibles from group

In [235]:
Image(filename= "math_model.png")
Out[235]:

Transmission model

Above an image is shown which summarizes the ODEs of [1] model in a scheme. The parameters in the model can be summarized as following:

  • $Sen$ and $Res$ indicate the antibitotic-sensitive and -resistand N. gonorrhoeae strains
  • $G = \{L,H\}$ is the set of low and high sexual activity groups, $i \in G$
  • Each sexual activity group $N_i$ consists of:
    • susceptible hosts $S_i$
    • hosts infected with an antibiotic-sensitive strain, $I_{Sen_i}$
    • hosts infected with an antibiotic-resistant strain, $I_{Res_i}$
  • Hosts are redistributed to either the same of the other sexual activity group at rate $\gamma$
  • Hosts leave or enter the population at rate $\alpha$
  • Susceptible hosts become infected depending on:
    • the partner change rate, $\pi_i$
    • the transmission probability per partnership, $\beta_{ij}$
    • the sexual mixing matrix $\rho_{ij}$
  • Hosts infected with an antibiotic-sensitive strain can:
    • recover spontaneously at rate $\nu$
    • receive treatment at rate $\tau$
  • Hosts receiving treatment recover at rate $\tau(1 - \mu)$
  • Hosts develop resistance during treatment with probability $\mu$
  • Host infected with an antibiotic-resistance strain can only recover spontaneously at rate $\nu$

The sexual mixing matrix describes how many partnerships are formed within and outside the host's activity group:

$$ \rho_{ij} = \epsilon \delta_{ij} + (1 - \epsilon) \frac{\pi_j N_j}{\sum_{k \in G} \pi_k N_k} $$

Here, $\delta_{ij} = 1$ if $i = j$ and 0 otherwise. Model parameters were estimated from sexual behavior data, calibrated trough model simulation or informed by literature. The parameters can be found on page 5 of [1].

The differential equations are solved below for two different groups: Men having sex with men (MSM) and heterosexual men (HMW).

In [236]:
## imports for the math model
import numpy as np
import matplotlib.pyplot as plt
import scipy as sp
from scipy.integrate import odeint
from math_model_functions import *

## imports for the network model
import networkx as nx
from networkx.utils.random_sequence import (powerlaw_sequence)
from numba import jit
from network_model import *

## imports for the stochastic network model
from stochastic_network_model import *
import matplotlib.lines as mlines

from matplotlib import animation, rc
from IPython.display import HTML
rc('animation', html='jshtml')
# import networkx.utils.create_degree_sequence as create_degree_sequence
%matplotlib notebook
# TODO: Change resolution of animation

MSM Parameters

Parameters for the "Men having sex with men group"

In [237]:
Nh = 5.3e-2
N = np.array([1 - Nh, Nh])

S = np.array([N[0]*0.995,N[1]*0.90])
Ires = np.array([0, 0])
Isen = np.array([N[0]*0.005, N[1]*0.1])

pi = np.array([0.41, 30.49])


BLL = np.array([0.42,0.59,0.77])
BHH = np.array([0.25,0.30,0.40])
alpha = 1 / 29


gamma = 1
epsilon = np.array([0.30,0.57,0.80])
phi = np.array([0.48,0.64,0.81])
mu = 1e-3
D = np.array([0.14,0.19,0.25])

y0 = np.concatenate((S, Ires, Isen))
# y0 /= np.sum(y0)
t = np.arange(0,30,0.001)
sol_msm = []
for i in range(3):
    beta = np.array([[BLL[i], 1], [1, BHH[i]]])
    beta[0, 1] = beta[1, 0] = np.sqrt(beta[0, 0] * beta[1, 1])
    sol_msm.append(odeint(math_model, y0, t, args=(N, pi, alpha, beta, gamma, epsilon[i], phi[i], mu, D[i])))
    
sol_msm = np.array(sol_msm)
In [238]:
# calculate R0 for low and high activity groups
R0_low = BLL[1] * pi[0] * D[1]
R0_high = BHH[1] * pi[1] * D[1]

# average R0 over population
R0 = Nh * R0_high + (1 - Nh) * R0_low
print(R0)
0.13563535699999998

HMW Parameters

Parameters for the heterosexual men

In [239]:
Nh = 6.3e-2
N = np.array([1 - Nh, Nh])

S = np.array([N[0]*0.995,N[1]*0.90])
Ires = np.array([0, 0])
Isen = np.array([N[0]*0.005, N[1]*0.1])

pi = np.array([0.25, 4.57])


alpha = 1 / 29
BLL = np.array([0.79,0.87,0.94])
BHH = np.array([0.63,0.72,0.81])

gamma = 1
epsilon = np.array([0.53,0.73,0.89])
phi = np.array([0.36,0.50,0.66])
mu = 1e-3
D = np.array([0.46,0.55,0.66])

y0 = np.concatenate((S, Ires, Isen))
# y0 /= np.sum(y0)
t = np.arange(0,30,0.001)
sol_hmw = []
for i in range(3):
    beta = np.array([[BLL[i], 1], [1, BHH[i]]])
    beta[0, 1] = beta[1, 0] = np.sqrt(beta[0, 0] * beta[1, 1])
    sol_hmw.append(odeint(math_model, y0, t, args=(N, pi, alpha, beta, gamma, epsilon[i], phi[i], mu, D[i])))
    
sol_hmw = np.array(sol_hmw)
In [240]:
# calculate R0 for low and high activity groups
R0_low = BLL[1] * pi[0] * D[1]
R0_high = BHH[1] * pi[1] * D[1]

# average R0 over population
R0 = Nh * R0_high + (1 - Nh) * R0_low
print(R0)
0.22610098500000003

Figures and plots

In [241]:
plt.figure()
# plt.plot(t, sol_msm[1,:, 0], label='$S_{low}$')
# plt.plot(t, sol_msm[1,:, 1], label='$S_{high}$')
plt.plot(t, sol_msm[1,:, 2], label='$I_{res, low}$')
plt.plot(t, sol_msm[1,:, 3], label='$I_{res, high}$')
plt.plot(t, sol_msm[1,:, 4], label='$I_{sen, low}$')
plt.plot(t, sol_msm[1,:, 5], label='$I_{sen, high}$')

plt.legend(loc='best')
plt.xlabel('t (years)')
plt.title("Epidemic curve MSM group")
plt.ylabel("Fraction of population")
# plt.ylim(0, 1)
plt.xlim(0, 10)
plt.grid()
plt.show()
In [242]:
plt.figure()
# plt.plot(t, sol_msm[1,:, 0], label='$S_{low}$')
# plt.plot(t, sol_msm[1,:, 1], label='$S_{high}$')
plt.plot(t, sol_hmw[1,:, 2], label='$I_{res, low}$')
plt.plot(t, sol_hmw[1,:, 3], label='$I_{res, high}$')
plt.plot(t, sol_hmw[1,:, 4], label='$I_{sen, low}$')
plt.plot(t, sol_hmw[1,:, 5], label='$I_{sen, high}$')

plt.legend(loc='best')
plt.xlabel('t (years)')
plt.title("Epidemic curve HMW group")
plt.ylabel("Fraction of population")
# plt.ylim(0, 1)
plt.xlim(0, 30)
plt.grid()
plt.show()
In [243]:
plt.figure()
ratio_sol_msm = (sol_msm[:,:, 2] + sol_msm[:,:, 3])/(sol_msm[:,:, 2] + sol_msm[:,:, 3] + sol_msm[:,:, 4] + sol_msm[:,:, 5])
plt.plot(t, ratio_sol_msm[1,:], label = "MSM", c = "g",linewidth=1)
plt.fill_between(t,ratio_sol_msm[0,:], ratio_sol_msm[2,:], alpha=0.5, facecolor = "g")
ratio_sol_hmw = (sol_hmw[:,:, 2] + sol_hmw[:,:, 3])/(sol_hmw[:,:, 2] + sol_hmw[:,:, 3] + sol_hmw[:,:, 4] + sol_hmw[:,:, 5])
plt.plot(t, ratio_sol_hmw[1,:], label = "HMW", c = "b",linewidth=1)
plt.fill_between(t,ratio_sol_hmw[0,:], ratio_sol_hmw[2,:], alpha=0.5, facecolor = "b")
plt.legend(loc='best')

plt.xlabel('t (years)')
plt.ylabel("Fraction of infected individuals")
plt.grid()
plt.ylim(0,1)
plt.show()

Simple Contact network

  • We used a simplified version of the rules in the mathematical model
  • The network is a random erdos renyi network with 500 nodes and an average degree of 10

The parameters

  • Spontaneous healing rate $\nu$
  • Treatment rate $\tau$
  • Hosts develop resistance during treatment with probability $\mu$
  • Transmision probability $\beta$

Build a network

In [264]:
# # create power law graph
# exp = 2.4
# seed = 123456789
# num = 1000
# graph = power_law_graph(exp, num, seed)

# or, create random network
degree = 10
num = 500
p = degree/num
seed = 1234
graph = random_network(p, num, seed)

# compute average degree
degrees = dict(graph.degree())
sum_of_edges = sum(degrees.values())
average_deg = sum_of_edges / len(degrees)
print(average_deg)

# compute r0 (TODO)
9.752
In [265]:
plt.figure()
pos = nx.spring_layout(graph)
nx.draw(graph, pos, node_color='b', node_size=10, with_labels=False)
plt.show()

Implementation of disease spreading

In [268]:
# set person class to each node
for i in range(len(graph)):
    graph.node[i]["Data"] = Individual(i)

# create set of infecteds, initialize with 10 infecteds
init = np.random.choice(len(graph), 10, replace=False)

# time steps
t = 0
t_tot = 20
time_step = 365
dt = 1 / (time_step)

# constants of model
D = np.array([0.19,0.55]) * time_step
mu = 1e-3
PHI = [0.65, 0.50]
BETA = [0.59, 0.87]

# keep track of infecteds and steps
num_steps = int(t_tot / dt)
num_infected = np.zeros((2, num_steps))
num_res = np.zeros((2, num_steps))


p_inf = 0.8
# iterate over sextypes (MSM and HMW)
for sextype in range(2):
    phi = PHI[sextype]
    beta = BETA[sextype]
    nu = (1 - phi)/D[sextype]
    tau = phi/D[sextype]
    beta = 1-(1-1/(average_deg*(1-p_inf)))**(nu)
    num_res[sextype,:], num_infected[sextype,:], graph = network_model(beta, tau, nu, mu, init, num_steps, copy.deepcopy(graph), True)
100%|██████████| 7300/7300 [00:44<00:00, 162.88it/s]
100%|██████████| 7300/7300 [00:42<00:00, 169.92it/s]

Plotting

In [280]:
# plot number of infecteds
t = np.linspace(0,t_tot, num_steps)    
plt.figure()
print("Mean number of infecteds:", np.mean(num_infected))
plt.plot(t, num_infected[0,:], label='all infecteds MSM')
plt.plot(t, num_res[0,:], label='resistant infecteds MSM')
plt.plot(t, num_infected[1,:], label='all infecteds HMW')
plt.xlabel("Time (years)")
plt.ylabel("number of people")
plt.plot(t, num_res[1,:], label='resistant infecteds HMW')
plt.legend()
plt.show()
Mean number of infecteds: 288.4821917808219
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-280-cb51f78d178f> in <module>()
      3 plt.figure()
      4 print("Mean number of infecteds:", np.mean(num_infected))
----> 5 plt.plot(t, num_infected[0,:], label='all infecteds MSM')
      6 plt.plot(t, num_res[0,:], label='resistant infecteds MSM')
      7 plt.plot(t, num_infected[1,:], label='all infecteds HMW')

IndexError: too many indices for array
In [281]:
plt.figure()
plt.plot(t,num_res[0,:]/num_infected[0,:], label='fraction of resistant strain MSM')
plt.plot(t,num_res[1,:]/num_infected[1,:], label='fraction of resistant strain HMW')
# plt.plot(num_res)
plt.xlabel("Time (years)")
plt.ylabel("Fraction of infected individuals")
plt.legend()
plt.show()
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-281-9a4093e5464e> in <module>()
      1 plt.figure()
----> 2 plt.plot(t,num_res[0,:]/num_infected[0,:], label='fraction of resistant strain MSM')
      3 plt.plot(t,num_res[1,:]/num_infected[1,:], label='fraction of resistant strain HMW')
      4 # plt.plot(num_res)
      5 plt.xlabel("Time (years)")

IndexError: too many indices for array

Animation

In [279]:
# set person class to each node
for i in range(len(graph)):
    graph.node[i]["Data"] = Individual(i)

# create set of infecteds, initialize with 10 infecteds
init = np.random.choice(len(graph), 10, replace=False)

# time
t = 0
t_tot = 5
time_step = 365
dt = 1 / (time_step)

# constants of model
D = np.array([0.19,0.55]) * time_step
mu = 1e-3
PHI = [0.65, 0.50]
BETA = [0.59, 0.87]

# choose MSM or HMW to model
sextype = 0
phi = PHI[sextype]
beta = BETA[sextype]
nu = (1 - phi)/D[sextype]
tau = phi/D[sextype]
beta = 1-(1-1/(average_deg*(1-p_inf)))**(nu)

# initial state
_, _, graph = network_model(beta, tau, nu, mu, init, 0, graph, doInit=True, disable_progress=True)

# get values of disease status
values = []
for j in range(len(graph.nodes())):
    if graph.node[j]["Data"].disease_status == 0:
        values.append("Black")
    elif graph.node[j]["Data"].disease_status == 1:
        values.append("Blue")
    elif graph.node[j]["Data"].disease_status == 2:
        values.append("Red")
        
# create figure and set axes
fig = plt.figure(figsize = (10,5))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
im1, = ax2.plot([], [], c = 'b', label='all infecteds MSM')
im2, = ax2.plot([], [], c = 'r', label='resistant infecteds MSM')
ax2.legend()
ax2.set_xlim(0,t_tot)
ax2.set_ylim(0,num)

# initialize
pos = nx.spring_layout(graph) 
plt.close()
ani_step = 365//5
frames = (t_tot*365)//ani_step
num_infected = np.zeros(ani_step*frames)
num_res = np.zeros(ani_step*frames)

# update function for animation
def update(i):
    if i<1:
        temp_res, temp_infected, _ = network_model(beta, tau, nu, mu, init, 1, graph, doInit=False, disable_progress=True)
    else:
        temp_res, temp_infected, _ = network_model(beta, tau, nu, mu, init, ani_step, graph, doInit=False, disable_progress=True)
    num_infected[i*ani_step:(i+1)*ani_step] = temp_infected
    num_res[i*ani_step:(i+1)*ani_step] = temp_res
    t = np.linspace(0,i*ani_step/365, i*ani_step) 
    im1.set_data(t, num_infected[:i*ani_step])
    im2.set_data(t, num_res[:i*ani_step])
    ax1.clear()
    
    # solution at step i
    values = []
    for j in range(len(graph.nodes())):
        if graph.node[j]["Data"].disease_status == 0:
            values.append("Black")
        elif graph.node[j]["Data"].disease_status == 1:
            values.append("Blue")
        elif graph.node[j]["Data"].disease_status == 2:
            values.append("Red")
    nx.draw(graph, pos=pos, node_color=values, node_size=10, ax=ax1, clim = 2)
    
    # Scale plot ax
    ax1.set_title("Frame %d:    "%(i+1), fontweight="bold")
    ax1.set_xticks([])
    ax1.set_yticks([])
    
# animate!
animation.FuncAnimation(fig, update, frames=frames, interval=100, repeat=True)
25
Out[279]:


Once Loop Reflect

Stochastic network model

The stochastic network model is based on an existing model by Kretzschmar, M. et al. (1996) [4]

The system consists of Persons defined by:

  • Gender (male / female)
  • Age (15 - 64)
  • Sexual activity (high / low)
  • Disease status (uninfected or (non-)resistant (a)symptomatic)

During each timestep, Partnerships are formed between two Persons. A Partnership is defined by:

  • Type (casual or steady)
  • The two Persons involved

The sexual contact network is therefore dynamic in time

Each timestep of the model consists of seven parts:

  1. Partnership formation
  2. Disease transmission
  3. Separation of partnerships
  4. Replacement
  5. Recovery
  6. Treatment
  7. Advance time

Initialize model

  1. Initialize 10,000 Persons of random age. 5% of Persons with age < 35 have high sexual activity
  2. Perform 2,000 iterations to obtain the steady state
  3. Infect 10% of the core group (high sexual activity) with non-resistant asymptomatic disease

The resistant disease enters the system through treatment; treatment will result in contracting the resistant disease with probability 0.01%

In [228]:
#initialize System
s = System()
s.set_constants()
s.initialize()

#initialize lists to store counts
singles_list = []
steady_list = []
partnerships_list = []

#start time steps to reach steady state
start = time.time()
for i in range(2000):
    s.time_step()
    singles_list.append(s.number_of_singles)
    steady_list.append(s.number_of_steady)
    partnerships_list.append(len(s.partnerships))

print("TIME TO INITIALIZE MODEL: %.2f"%(time.time() - start))

#print core size (sexual_activity == 1)
core_size = len(s.high_sexual_activity)
print("Core size: %d"%(core_size))

#random infection of the core group
infect = np.random.choice(list(s.high_sexual_activity),int(0.1 * core_size))
for i in infect:
    p = s.persons[i]
    p.asymptomatic_infection(s)

#set System time to zero
s.time = 0
TIME TO INITIALIZE MODEL: 53.32
Core size: 220

Time-evolution towards the steady state

In [229]:
# #print percentage of singles in the System
# singles = len(s.singles)/len(s.persons) * 100
# print("Percentage of singles: %f"%(singles))

# #print percentage having one partner
# one_partner = 0
# for p in s.persons:
#     if len(p.current_partners) == 1:
#         one_partner += 1
# one = one_partner/len(s.persons) * 100
# print("Percentage having one partner: %f"%(one))
# print("Percentage having more than one partner: %f"%(100 - singles - one))
In [230]:
%%capture
#plot change over time, reaching the steady state
fig, ax = plt.subplots(figsize=(9,6))
ax.plot(singles_list,label="number of singles")
ax.plot(steady_list,label="number of steady partnerships")
ax.plot(partnerships_list,label="number of partnerships")
ax.plot(len(s.persons)-np.array(singles_list),label="number of persons in a partnership")
ax.plot(np.array(partnerships_list) - np.array(steady_list),label="number of casual partnerships")
ax.set_title("Time-evolution of Persons and Partnerships during initialization")
ax.set_xlabel("Days")
ax.set_ylabel("Counts")
ax.legend();
In [231]:
fig

Evolve system over time

In [232]:
#initialize lists
singles_list = [s.number_of_singles]
steady_list = [s.number_of_steady]
partnerships_list = [len(s.partnerships)]
symptomatic_list = [s.number_of_symptomatic]
asymptomatic_list = [s.number_of_asymptomatic]
res_symp_list = [s.number_of_res_symp]
res_asymp_list = [s.number_of_res_asymp]

#store starting time
start = time.time()

for i in range(12000):
    #make a single time step
    s.time_step()
    
    #append new values to lists
    singles_list.append(s.number_of_singles)
    steady_list.append(s.number_of_steady)
    partnerships_list.append(len(s.partnerships))
    symptomatic_list.append(s.number_of_symptomatic)
    asymptomatic_list.append(s.number_of_asymptomatic)
    res_symp_list.append(s.number_of_res_symp)
    res_asymp_list.append(s.number_of_res_asymp)

#print evaluation time
print(time.time() - start)
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-232-edf646322973> in <module>()
     13 for i in range(12000):
     14     #make a single time step
---> 15     s.time_step()
     16 
     17     #append new values to lists

~/Documents/GitHub/CSS2018/Code/stochastic_network_model.py in time_step(self)
    693             if p.gender == 0 and p.time_since_infection > self.incubation_time_men:
    694                 p.symptomatic_infection(self)
--> 695             elif p.gender == 1 and p.time_since_infection > self.incubation_time_women:
    696                 p.symptomatic_infection(self)
    697 

KeyboardInterrupt: 

Time-evolution of the population

In [ ]:
%%capture
#show change in counts over time
fig, ax = plt.subplots(figsize=(9,6))
ax.plot(singles_list,label="number of singles")
ax.plot(steady_list,label="number of steady partnerships")
ax.plot(partnerships_list,label="number of partnerships")
ax.plot(len(s.persons)-np.array(singles_list),label="number of persons in a partnership")
ax.plot(np.array(partnerships_list) - np.array(steady_list),label="number of casual partnerships")
ax.set_title("Evolution of the system over time")
ax.set_xlabel("Time in days")
ax.set_ylabel("Counts")
ax.legend();
In [ ]:
fig

Time-evolution of infecteds over time

In [ ]:
%%capture
#show change in infecteds over time
fig,ax = plt.subplots(figsize=(9,6))
ax.plot(symptomatic_list,label="symptomatically infected")
ax.plot(asymptomatic_list,label="asymptomatically infected")
ax.plot(res_symp_list,label="resistant symptomatically infected")
ax.plot(res_asymp_list,label="resistant asymptomatically infected")
ax.set_title("Evolution of infecteds over time")
ax.set_xlabel("Time in days")
ax.set_ylabel("Counts")
ax.legend();
In [ ]:
fig

Fraction of resistant infecteds over all infecteds over time

In [ ]:
%%capture
#show change in resistant infections over time
res_inf = np.array(res_symp_list) + np.array(res_asymp_list)
inf = res_inf + np.array(symptomatic_list) + np.array(asymptomatic_list)
fig,ax = plt.subplots(figsize=(9,6))
ax.plot(res_inf/inf)
ax.set_title("Fraction of resistant infecteds over all infecteds over time")
ax.set_xlabel("Time in days")
ax.set_ylabel("Fraction of resistant infections");
In [ ]:
fig

Infectiousness of the disease ($R_0$)

Although the disease might seem quite infectious, simulations have shown $R_0=0.90 \pm 2.93$ (N=1,000, 95% CI). The spread of the disease is very dependent on where the disease is started. If a person with high sexual activity is infected, it will quickly spread throughout the entire core group whereas a person with low sexual activity may only affect a select few different partners

Animation

In [210]:
%%capture
#initialize System
s = System()
s.set_constants()
s.initialize()

#evolve to steady state
for i in range(2000):
    s.time_step()
    
#random infection of the core group
core_size = len(s.high_sexual_activity)
infect = np.random.choice(list(s.high_sexual_activity),int(0.1 * core_size))
for i in infect:
    p = s.persons[i]
    p.asymptomatic_infection(s)

#dictionary mapping disease_status to a color
color_dict = {0:'black',1:'blue',2:'green',3:'red',4:'orange'}
dis_dict = {0:'Uninfected',1:'Non-resistant symptomatic',
            2:'Non-resistant asymptomatic',3:'Resistant symptomatic',
            4:'Resistant asymptomatic'}

def plot_graph(i):
    """
    Draw graph for a single time step
    """
    global s, pos, ax
    
    if i > 0:
        for j in range(365):
            s.time_step()
    
    persons = set()
    edges_list = []
    for identifier in s.high_sexual_activity:
        p = s.persons[identifier]
        for partner in p.current_partners:
            persons.add(partner)
            persons.add(identifier)
            edges_list.append((identifier,partner))
    persons_list = sorted(list(persons))

    sizes = []
    disease = []
    for identifier in persons_list:
        disease.append(color_dict[s.persons[identifier].disease_status])
        sizes.append(15+30*s.persons[identifier].sexual_activity)

    ax.clear()
    G = nx.Graph()
    G.add_nodes_from(persons_list)
    G.add_edges_from(edges_list)
    nx.draw(G,node_color=disease,node_size=sizes,ax=ax)
    handles = []
    for index in range(5):
        handles.append(mlines.Line2D([],[],color=color_dict[index],marker='o',linestyle='None',
                                    markersize=10,label=dis_dict[index]))
    ax.legend(handles=handles,loc='upper left')
    ax.set_title("Core group and their partners\nT = %d years"%i)

fig,ax = plt.subplots(figsize=(10,5))
anim = animation.FuncAnimation(fig,plot_graph,frames=34,interval=1000,repeat=False)
In [211]:
anim
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-211-027b5d0310db> in <module>()
----> 1 anim

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/IPython/core/displayhook.py in __call__(self, result)
    255             self.start_displayhook()
    256             self.write_output_prompt()
--> 257             format_dict, md_dict = self.compute_format_data(result)
    258             self.update_user_ns(result)
    259             self.fill_exec_result(result)

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/IPython/core/displayhook.py in compute_format_data(self, result)
    149 
    150         """
--> 151         return self.shell.display_formatter.format(result)
    152 
    153     # This can be set to True by the write_output_prompt method in a subclass

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/IPython/core/formatters.py in format(self, obj, include, exclude)
    178             md = None
    179             try:
--> 180                 data = formatter(obj)
    181             except:
    182                 # FIXME: log the exception

<decorator-gen-9> in __call__(self, obj)

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/IPython/core/formatters.py in catch_format_error(method, self, *args, **kwargs)
    222     """show traceback on failed format call"""
    223     try:
--> 224         r = method(self, *args, **kwargs)
    225     except NotImplementedError:
    226         # don't warn on NotImplementedErrors

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/IPython/core/formatters.py in __call__(self, obj)
    343             method = get_real_method(obj, self.print_method)
    344             if method is not None:
--> 345                 return method()
    346             return None
    347         else:

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/matplotlib/animation.py in _repr_html_(self)
   1484             return self.to_html5_video()
   1485         elif fmt == 'jshtml':
-> 1486             return self.to_jshtml()
   1487 
   1488 

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/matplotlib/animation.py in to_jshtml(self, fps, embed_frames, default_mode)
   1467                 self.save(f.name, writer=HTMLWriter(fps=fps,
   1468                                                     embed_frames=embed_frames,
-> 1469                                                     default_mode=default_mode))
   1470             # Re-open and get content
   1471             with open(f.name) as fobj:

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/matplotlib/animation.py in save(self, filename, writer, fps, dpi, codec, bitrate, extra_args, metadata, extra_anim, savefig_kwargs)
   1257                     for anim, d in zip(all_anim, data):
   1258                         # TODO: See if turning off blit is really necessary
-> 1259                         anim._draw_next_frame(d, blit=False)
   1260                     writer.grab_frame(**savefig_kwargs)
   1261 

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/matplotlib/animation.py in _draw_next_frame(self, framedata, blit)
   1294         # post- draw, as well as the drawing of the frame itself.
   1295         self._pre_draw(framedata, blit)
-> 1296         self._draw_frame(framedata)
   1297         self._post_draw(framedata, blit)
   1298 

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/matplotlib/animation.py in _draw_frame(self, framedata)
   1812         # Call the func with framedata and args. If blitting is desired,
   1813         # func needs to return a sequence of any artists that were modified.
-> 1814         self._drawn_artists = self._func(framedata, *self._args)
   1815         if self._blit:
   1816             if self._drawn_artists is None:

<ipython-input-210-e04a745e332f> in plot_graph(i)
     29     if i > 0:
     30         for j in range(365):
---> 31             s.time_step()
     32 
     33     persons = set()

~/Documents/GitHub/CSS2018/Code/stochastic_network_model.py in time_step(self)
    645                 if np.random.random() < self.recovery_rate_asymptomatic_men:
    646                     p.cure(self)
--> 647             elif p.gender == 1 and p.disease_status in [1,3]:
    648                 if np.random.random() < self.recovery_rate_symptomatic_women:
    649                     p.cure(self)

KeyboardInterrupt: 

Determine $R_0$

In [ ]:
#import to perform in parallel
from joblib import Parallel, delayed

#Note that R0 is very high, due to treatment. If Person A is asymptomatic and has a Partnership with Person B,
#who contracts the symptomatic disease, Person B gets treatment. Person A is however still infectious, so Person
#B will almost immediately contract the disease again. To 'solve' this, we convert the list of Person identifiers
#to a set, such that infections are counted uniquely per Person.

def get_R0():
    """
    Gets a single sample of R0 from the System    
    """
    #initialize System
    s = System()
    s.set_constants(r0=True)
    s.initialize()

    #evolve to steady state
    for j in range(1000):
        s.time_step()
    
    #infect Person 0
    s.persons[0].asymptomatic_infection(s)
    
    #loop until Person 0 is cured
    while not s.persons[0].disease_status == 0:
        s.time_step()
    
    #return R0
    return s.r0_infected

#perform 4 jobs in parallel
a = Parallel(n_jobs=4,verbose=5)(delayed(get_R0)() for i in range(1000))

#convert lists to sets, such that we only have the unique Persons infected
a = [len(set(x)) for x in a]

#print R0 with 95% CI
print(a)
print("R_0 = %f +/- %f"%(np.mean(a),1.96*np.std(a)))

References

[1] Fingerhuth, S. M., Bonhoeffer, S., Low, N., & Althaus, C. L. (2016). Antibiotic-resistant Neisseria gonorrhoeae spread faster with more treatment, not more sexual partners. PLoS pathogens, 12(5), e1005611.

[2] De, P., Singh, A. E., Wong, T., Yacoub, W., & Jolly, A. M. (2004). Sexual network analysis of a gonorrhoea outbreak. Sexually transmitted infections, 80(4), 280-285.

[3] Bansal, S., Grenfell, B. T., & Meyers, L. A. (2007). When individual behaviour matters: homogeneous and network models in epidemiology. Journal of the Royal Society Interface, 4(16), 879-891.

[4] Kretzschmar, M., van Duynhoven, Y. T., & Severijnen, A. J. (1996). Modeling prevention strategies for gonorrhea and chlamydia using stochastic network simulations. American Journal of Epidemiology, 144(3), 306-317.